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> ■ Abstract 
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A standard drift-diffusion model of space charge wave propagation in semicon- 
^ , ductors has been studied numericahy and analytically under dc voltage bias. 

Q \ For sufficiently long samples, appropriate contact resistivity and applied volt- 

age - such that the sample is biased in a regime of negative differential resis- 
tance - we find chaos in the propagation of nonlinear fronts (charge monopoles 
of alternating sign) of electric field. The chaos is always low-dimensional, but 
\ has a complex spatial structure; this behavior can be interpreted using a finite 

Q \ dimensional asymptotic model in which the front (charge monopole) positions 

J-^ ' and the electrical current are the only dynamical variables. 
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I. INTRODUCTION 



The dynamics of propagating solitary waves (pulses) and monotone fronts in nonlinear 
partial differential equations {pde^s) have been the focus of extensive research. Such problems 
are of interest in a wide range of fields including biology (population dynamics) , chemical 
reactions and combustion 0,0], plasma physics 0], and semiconductor electronic transport 
A common problem is to understand and predict the form and speed of the waves, as 
well as the possibility of multiple excitations in a sample of finite or infinite extent, and their 
interactions. We focus here on a model of electrical conduction in extrinsic semiconductors 
(involving time and only one spatial dimension) which exhibits negative differential resis- 
tance (NDR) and moving domains of high electric field. The model is specifically relevant 
to experiments on cooled bulk p-type Ge under voltage bias conditions, but much of 
the observed qualitative behavior applies to a broad class of semiconductor systems with 
space charge instabilities. Phenomena observed for the p-Ge system include time-periodic 
oscillation of the current in a purely resistive external circuit under dc voltage bias due 
to the periodic creation of a solitary wave at the injecting contact, its motion inside the 
semiconductor and its annihilation at the receiving contact [^. There is some similarity to 
the Gunn effect in n-GaAs, except that: (i) the local current density versus field characte- 
ristics [see j{E) below] in p-Ge presents an increasing third branch after the NDR region, 
and (ii) the solitary waves in p-Ge move much more slowly than the carrier drift velocity 
(the case in the usual Gunn effect) due to the generation-recombination dynamics of ionized 
traps which dominate the transport properties. Other experimental observations include 
intermittency near the onset of the oscillatory instability PJ^, and "spatiotemporal" chaos 
under combined dc and ac voltage bias [1^. Another important feature is the integral con- 
straint which corresponds to voltage bias applied across sample contacts in semiconductor 
problems. Integral constraints also occur in other situations, for example, expressing mass 
conservation in problems of phase separation in binary mixtures and in certain biological 
problems 111,0,0. 

Many of these phenomena have been successfully explained by means of a drift-diffusion 
model which includes impurity trapping of mobile holes and impact ionization of neutral 
acceptors |]6|,p!3|. Although much work has been done on this model problem (see and 



references therein), important basic questions concerning its asymptotic description and 
chaos under time independent voltage bias are still open. In this paper, we present numerical 
simulation results which show chaos under dc voltage bias associated with multiple shedding 
of wavefronts. Multiple shedding of wavefronts occurs for appropriate values of contact 



resistivity, and was recently predicted on the basis of asymptotic calculations []15[. Here we 
introduce a finite dimensional model which provides a simplified description of space charge 
wave dynamics in long samples. This model uses relevant information from the asymptotics 
in Ref. although we do not rigorously derive it from such asymptotic calculations. 
Nevertheless, solutions of the simplified model are in good agreement with the results of 
direct numerical simulations. 

The paper is organized as follows. In Section II we present a drift-diffusion model which 
accurately describes the propagation of space charge wavefronts in extrinsic semiconductors 
under dc voltage bias such as p-Ge. In Section III we present numerical simulation results 
of the drift- diffusion model which indicate that wavefronts may propagate chaotically due to 
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multiple shedding of wavefronts for appropriate bias in sufficiently long samples. Section IV 
presents the asymptotic model and numerical simulations thereof which are then compared 
with results from the drift-diffusion model. Conclusions are finally presented in Section V. 



II. THE REDUCED DRIFT-DIFFUSION MODEL 



In dimensionless form, the drift-diffusion model equations for a sample of length L can 
be written in a form p!6| , p^ : 



d'^E K + R / V dE 
dxdi^ \K + R^ 

1 

E{x,t) dx = (p, (2) 



L 







E{0,t) = pJ{t). (3) 



The first equation describes the spatio-temporal evolution of the electric field E{x,t) 
inside the sample, where J{t) is the total current density. The transport coefficients V, K 
and R are, respectively, the average (drift) velocity of microscopic charge carriers (holes in 
the case of p-Ge), and coefficients describing the creation of free carriers via impurity impact 
ionization and the destruction of free carriers by capture onto an available impurity trapping 
site (i.e., neutral acceptor). All are nonlinear functions of electric field, and their forms and 
plots have been discussed extensively in the literature In this paper we use the 



same forms as in Ref. [Q. Equation (|^) is a global constraint which expresses the voltage 
bias condition, and Eq. (|]) is a boundary condition which represents the Ohmic injecting 
contact at X = with contact resistivity p > 0. We refer to Eqs. (|l|) - (||) as the reduced drift- 
diffusion model because they are derived from a full drift- diffusion model by systematically 
dropping terms that correspond to short length and time scale processes of diffusion and 
displacement current, respectively. For a precise derivation as well as a complete table of 
conversion factors to dimensional units, see |^ and |T^. Some of them are: time 2.1 10" 



ms, length 0.01 mm, electric field 10 V/cm, density current 128.16 mA/cm^, cross-sectional 
area 0.16 cm^. We have adopted the same symbol for both dimensional and non-dimensional 
variables. 

The qualitative nature of much of the dynamical behavior found in the reduced drift- 
diffusion model — e.g., the instability of the stationary electric-field profile and propagating 



high-field domains |[T7|-[T9| — depends only on the presence of a region of negative slope of 



the homogeneous stationary current density i{E) = V {E) {aK{E) /[K[E) + R{E)\ — 1}, 



1^ , P^J18| over an interval of positive fields, and not on the exact form of the underlying 
coefficients. This is particularly true when the sample is closely compensated (the ratio 
of the acceptor concentration to the donor concentration, a, is only slightly larger than 
1). Then j{E) is N-shaped for large enough positive fields: there is an interval [Em, Em) 
between the abcissas of the maximum [J^Em) = Jm] and the minimum [j{Em) = jm > 0] 
current density for which dj/dE < and j{E) > 0, as shown in Fig. |l|. Also shown in 
Fig. ID is the injecting contact characteristic which plays a crucial role in determining when 
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new fronts are injected into the sample. The critical current density, Jc, corresponds to the 
intersection of the contact characteristic with the homogeneous stationary current density. 
The role of Jc has been elucidated using a rigorous asymptotic analysis of the system Eqs. 
(HI) - (3) in the limit as L — > oo |jl5|, and is discussed further in Section IV. 



III. NUMERICAL RESULTS 

To solve the system of equations (||) - (0) for E and J, we discretize the equations 
using finite difference approximations to the derivatives and employ an implicit method to 
generate the solution. The initial condition for the electric field is spatially uniform with 
a value that is consistent with the global constraint, Eq. @. In Fig. ^ we show a space- 
time plot of the electric field and associated current density for cf) = 6.25 V/cm, just above 
the threshold voltage value for which propagating domain behavior occurs. The gray scale 
ranges from 5.3 V/cm (black) to 14.2 V/cm (white), and this scale is used in all similar plots 
that follow. The dimensionless sample length is 3800 corresponding to a real p-Ge sample 
of length 3.87 cm, and the contact resistivity p is 780 Qcm corresponding to a value of 
10.0 in dimensionless units. This case corresponds well to experimental data, but published 
data in p-Ge were only presented for one sample and relatively low bias values and contact 
resistivity. We clearly see that a single domain moves across the sample at constant speed, 
until it reaches the receiving contact. As it disappears, a new wave is created at the injecting 
contact and the process repeats periodically. The current versus time plot indicates that the 
current is steady when the domain moves in the sample interior, while there is an increase 
when the domain reaches the receiving contact. It is important to note that the fronts of 
changing electric field (or equivalently, the regions of nonzero charge density) are sharp in 
space relative to other physical length scales for this problem, i.e, the extent of the flat top 
domains or the sample length. We have found that this separation of length scales increases 
with sample length and holds for most biases of interest. It is only for voltages near the 
onset point that one tends to observe rounded solitary waves rather than well-separated 
pairs of fronts; this is the dominant space charge wave structure observed in shorter samples 
and has been extensively reported [0,|1^,|1^ . 



As the bias increases, the propagating domain becomes fatter and eventually a second 
small domain is nucleated and propagates part way into the sample; but it dies before 
reaching the receiving contact or merging with the larger domain. At even larger bias values 
the second domain merges with the primary domain near the receiving contact and this 
situation is shown in Fig. |^ which corresponds to = 7.25 V/cm. Again the current is 
plotted on the far right of the figure. When the first domain reaches the receiving contact 
the current increases. Instead of immediately starting the nucleation of a new wave, the area 
lost by the dying wave is gained by the trailing wave - note that the width of the trailing 
wave increases after the leading wave starts to disappear. The current increases, reaching 
a local maximum just before the trailing domain touches the leading domain, that is, the 
fronts collide. The current increases abruptly after the front collision and it rises to a global 
maximum at which point a new domain begins to nucleate at the injecting contact. As the 
domain forms, the current decreases and reaches a minimum at which point a new domain 
detaches and begins to propagate. Then the current increases until a second smaller domain 
is nucleated. Finally, the current settles to a rather low constant level as the two domains 



4 



move steadily and in unison across the interior region of the sample. Current behavior is 
apparently dominated by the major events involving the fronts: collisions with the contacts 
or with each other. This suggests the viability of a dynamical model that focuses on discrete 
front motions and the current J{t). 

At larger biases, the portion of the sample occupied by the high field value is larger, 
reducing the separation between domains, and giving more complicated E{x,t) structure 
and J — t behavior. In Fig. ^ we show a space-time plot and current for what appears to 
be a chaotic state for an applied bias of = 10.0 V/cm. The spatiotemporal dynamics 
possess a great deal of structure and complexity. The process of multiple domain shedding 
is similar to that for the previous case. The large current peaks correspond to nucleation 
of leading domains. The leading domains cross the sample without catching up or forward- 
colliding with any other high field regions, and are indicated by dark regions that extend all 
the way across the space-time diagram. Note also the larger spatial extent of the leading 
domains than in Figs. ^ and ^ The aperiodicity of the current is reflected in the irregular 
appearance of the maximum current peaks or, equivalently, of the dark strips that extend 
across the entire sample. In between them are a number of local maxima corresponding to 
the shedding of trailing domains. 

Figure ^ (a) shows a bifurcation diagram in which we plot all values of successive current 
maxima as a function of (p. An important feature in this diagram is the apparent presence 
of windows of chaotic behavior with a large number of points being visited. For periodic 
states we see a small number of points corresponding to perfectly repeating current maxima. 
Also, in the periodic regimes there are points where various branches merge or disappear and 
these correspond to the development or destruction of trailing domains. For the parameter 
values selected here we do not observe period doubling. We conjecture that the route to 
chaos here is of boundary crisis type in which the attractor collides with a periodic orbit on 
its basin boundary pO| . 

In Fig. 1^ (b) we show the largest Lyapunov exponent Ai versus (j) for the reduced drift- 
diffusion model. This unambiguously confirms the presence of chaos in the "chaotic" win- 
dows. The next two exponents have been calculated and are never positive, so that the chaos 
we see is of a low- dimensional variety. To compute the exponents, we used an algorithm 
outlined in Ref. ||2l| adapted for use with partial differential equations and using adaptive 
control of the integration time step The values of Ai are zero in the periodic regimes 
as they should be for periodic behavior. The smallness of Ai in the chaotic regime is easily 
understood by recalling the period of the system, about 1000 non-dimensional time units. 
This indicates that the chaos originates in processes that occur on time scales on the order 
of the front transit time across the sample. To our knowledge, this is the first time that 
chaos due to multiple shedding of wavefronts has been observed in a drift-diffusion model of 
this type. This type of chaotic behavior has not been reported for experiments on p-Ge with 
time-independent voltage bias, most likely because experimentally studied samples were too 
short, biases were not sufficiently large, or contact resistivity was too low. However, in 
early experimental studies of the Gunn effect in GaAs, J. B. Gunn observed that for 
long samples current oscillations were almost completely random, resembling white noise. 
He also found that short samples produce aperiodic oscillations when circuit impedance is 
sufficiently large. It is plausible that Gunn may have observed a similar form of chaos to 
the one we have found numerically. In the next section we develop an asymptotic model in 
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which the chaotic dynamics is understood to arise from the aperiodic nucleation of fronts at 
the injecting contact. 



IV. ASYMPTOTIC MODEL 



The asymptotic model used in this paper consists of describing the evolution of the cur- 
rent when all wavefronts are detached from the injecting contact by an appropriate ordinary 
differential equation (ode) for J, tracking the position of the wavefronts and proposing a 
simplified mechanism for creation and destruction of wavefronts. While our new model is 
compatible with the asymptotic calculations of Ref . |]15| , we have not rigorously derived it 
from these calculations. Instead, we have proposed a simplified dynamics to account for our 
numerical observations motivated by asymptotic results. 

We begin by assuming that Em < (p < in Fig. |I] and that Em/ Jm < P < E^/ jm- 
Then a Gunn effect mediated by solitary waves occurs as shown by the numerical 
simulation in Fig. |[ For appropriate parameter values, there appears a regular oscillation 
of the current caused by repeated creation, motion and destruction of high-field domains 
in the sample. High-field domains are formed by two wavefronts separating a region where 
the electric field is uniform and large from regions of uniform low field. Clearly, there are 
positively and negatively charged wavefronts, having dE/dx > or dE/dx < 0, respectively. 
Near the contacts, there are narrow boundary layer regions where the electric field changes 
abruptly. Creation of high-field domains occurs at the injecting contact, via an instability 
of the boundary layer which expels a high-field domain from the injecting contact to the 
bulk of the sample. Typically, the total current changes most during wavefront creation 
and destruction events. In the limit as L ^ oo, space and time scales are x/L and t/L, 
respectively |T^. Then j{E) = J, except in wavefronts and boundary layers at the contacts. 



If the field profile consists of a single high-field domain detached from the contacts, we 
have E = E^lJ) inside the domain and E = Ei{J), outside, where Ei < E2 < E^ are 
the three zeros of j{E) — J for jm < J < 3m, High and low-field regions are joined 

by wavefronts which are the unique solution of Eq. ([l|) (with zero right hand side) in the 
moving coordinate x = ^ ~ X±if)i dX±/dt = c±{J) (the signs -|- or - refer to the charge 
inside the wavefront) and appropriate boundary conditions. For example, at a positively 
charged wavefront, E Ei{J) as x ^ — cxo and E —>■ E^[J) as x — ^ +00. The numerically 
determined values of c±{J) are shown in Fig. ^. 

Boundary layers obey (most of times) a quasistationary version of Eq. ([^) with appropriate 
boundary conditions on a semiinfinite spatial support. The instantaneous value of the cur- 
rent J {t/L) determines the field profile in the low and high uniform- field regions and the 
velocity of the wavefronts. 

Next, assume that we have an initial field profile consisting of high-field domains 
(solitary waves), each formed by two wavefronts located at XY{t) < X^\t). We shall 
number the wavefronts so that X^^ (t) > X^'^^^ (t) , and if necessary we shall consider xi^"* = 
L and X^^ = 0. Then the positions X^\t) are given by 

X^\t)= r cM{s))ds, (4) 
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where t^^\_ denotes the time at which the i-th monopole (with positive or negative charge) 
was born at x = 0. 

The evolution of the total current density is determined by the bias condition (Q), which 
may be approximately evaluated as 

<P = E,{J) + [E,{J) - E,{J)] Y: - ^ (5) 

1=1 

(terms of order 1/L and smaller have been ignored here; note that the x'^ / L are of order 
unity). We can get an ode for J by differentiating Eq. (|^) with respect to time and then 
substituting dX^ /dt = c±{J) in the result. We obtain 

dJ 1 [E^-E,] 



dt L ^ + ^ 



n+c+ — n_cJ) (6) 



dxt^ , dX^^ 



where i goes from 1 to A^. The quantities n+ and n_ are, respectively, the number of positive 
and negative monopoles detached from the contacts (i.e., excluding possible monopoles at 
a; = and x = L), while j[ and jg denote the derivative of the static j{E) characteristic 
with respect to electric field, evaluated at Ei and E^, respectively. Notice that the system 
of equations (§) to (0) completely specifies the behavior of current and field profile on the 
scales x/L and t/L, except that we do not have conditions for determining when new fronts 
are emitted from the injecting contact. 

We start with the simple case of Fig. ^ the motion of a single high-field domain far 
from the contacts. J satisfies Eq. (^ with n+ = n_ = 1, i.e., dJ/ds = A{J) [c+(J) — c_(J)], 
where s = t/L and A{J) > 0. As shown in Figure]^, c+(J) [resp. c_(J)] is a decreasing 
(resp. increasing) function of the current. Therefore, J evolves exponentially fast toward 
the zero of the right hand side of this equation, J = J*. When the leading wavefront at 
X = X_ arrives at x = L, it disappears almost instantaneously in the scale s, and we 
obtain dJ/ds = A{J) c^{J) > 0, so that the current increases. The injecting boundary layer 
near x = ceases to be quasistationary when J surpasses the value Jc at which the line 
J = E/p intersects the second (decreasing) branch of J = j{E). The precise description 
of the instability, which results in expelling a narrow high-field domain from a; = to the 
interior of the sample, can be found in ||15[. It is enough to say that a certain semiinfinite 
problem has to be solved numerically and matched to the resulting situation with a narrow 
high-field domain (consisting of a region of = E^{J) bounded by positively and negatively 
charged wavefronts) near x = and a high-field region from a positively charged wavefront 
to X = L. In the new situation, we have dJ/ds = A{J) [2c+(J) — c_(J)], and J tries to go 
toward the zero J = J+ of 2c+(J) — c_(J). Depending on the resistivity p, Jc can be larger 
than 7+ (and then J decreases toward J+), or e (J*, J+) (and then J starts increasing, 
and a second high-field domain may be expelled from x = 0). The simplest case, Jc > J~^, 
was described asymptotically in [T^. Provided cj) is large enough, J evolves exponentially 
fast towards J^. When the old domain leaves the sample, only two wavefronts (bounding 
the new high-field domain) remain, and J evolves exponentially fast toward J*, so that a 
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period of the oscillation is completed; see Fig. 0. The second case, Jc G (J*, J"*") is more 
complicated: numerical simulations show that multiple high-field domains may coexist in 
the sample at the same time as in Fig. ^. 

To achieve a simplified description of the current oscillation, valid for any positive value 
of Jci we proceed to examine further the shedding process. The simplest rule to determine 
when a new wavefront is shed from the injecting contact would be as follows: a positive 
(negative) front is emitted at the instant that J passes through with positive (negative) 
time derivative. However, this rule neglects the time needed for sufficient charge to be 
injected at the contact to form a propagating front. We may estimate the effective delay 
time by considering the time evolution of Eq. (|I]) evaluated at the injecting contact, 

. ,K + R 
u + J — — — u = 

+ J^^\J-3{E)i (8) 

where u{t) = dE{0, t)/dx and the argument of V, K and R is i?(0, t) = pJit), i.e., the value 
of electric field at the injecting contact. In this equation, we can think of J{t) as driving 
charge injection processes which determine front launching. Based on extensive simulations 
of the reduced model, we have found that u{t) must attain a sufficiently large positive or 
negative value for the front to detach and begin to propagate. This value is generally found 
to lie between 50% and 90% of the steady state value that \u\ would have in the case of 
no propagating fronts, i.e., where E near the injecting contact rapidly rises to the Es value 
or rapidly falls to the Ei value. The asymptotic model system is fully defined once the 
threshold is set and consists in Eqs. (^) to (|]). Thus, in the limit of an infinitely long 
sample, terms of order e = 1/L drop and we arrive at a low-dimensional dynamical system, 
which consists essentially of: (i) propagating negative and positive charge points that move 
according to the Eq. (H), (ii) subject to the conservation law Eq. (|^), (iii) which produce a 
measurable current according to Eq. (|), and are created according to Eq. (H). Note that 
in the e — limit the current will exhibit slope discontinuties at the formation times t^\_ 

and destruction times tl^\ but will be otherwise continuous and governed by Eq. (^). We 
note that the rigorous foundations of this and similar asymptotic models have been explored 



recently using singular perturbation methods |r5|,|25 



In this paper, we estimate the order-of-magnitude of the time delay associated with 
wavefront formation by evaluating Eq. (|^) for J J^. This implies a relaxation time of 

(9) 



^ - (a - l)KV' 



i.e., the approximate time for u{t) to go from a value of m(0) = —1 to u{t) = 0. By dJc/dt, 
we refer to the value of dJ/dt when J crosses Jc- Then, we adopt the criterion that a 
new front is born at x = at the time t + ar where t is the time at which J = Jc, and 
ar is a delay time. Here a is a number of order one which is determined from simulation 
of the reduced model for a particular bias voltage and then assumed to apply over the 
complete range. Wavefront destruction is assumed occur instantaneously at times t^*^ when 

ignore the finite duration of (fast) monopole 
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destruction stages which is equivalent to the well-justified approximation of neglecting the 



diffusive boundary layer at the receiving contact ||T5[. It should also be kept in mind that 
the index instantaneously decreases by one when wavefronts downstream collide with one 
another. 

Let us now use the above asymptotic model to interpret the simulation results for Eqs. 
(S) " dl)- The case of contact resistivity such that Jc > J"*" has been explained already: we 
obtain the usual Gunn effect with at most one solitary wave detached from the contacts for 
any time |lT5| , p5| ; see Fig. |. Let us assume now that the contact resistivity is such that 
Jc G (J*, J"*"). Then the current will increase after creation of a solitary wave, because 
2c+ — c_ > and multiple wave shedding is possible [jT5|. This situation is shown in Fig. 
1^, which shows simulation results of our simplified asymptotic model for similar parameters 
to those of Figure ^. The latter is depicted using data from direct numerical simulation of 
the reduced pde model Eqs. (|l]) - (|). To obtain Fig. ^ the values of a were set to 13.75 for 
positive front emission and 7.82 for negative front emission. We use these same values in 
the data of Fig. ||, which has same bias values as Fig. |[ The chaos appears to be closely 
tied to the asynchronous emission of fronts. This explains why the chaos observed for this 
partial differential equation system is low- dimensional. It is interesting to speculate how the 
maximal number of domains possible might scale with system size and contact resistivity. 



V. CONCLUSION 



We have utilized asymptotic analysis of a pde model (which describes the trap-dominated 
slow Gunn effect in a long sample) to explain the dynamics of space charge waves and cur- 
rent vs. time, including low-dimensional chaos which is nonetheless accompanied by spatially 
complex structure suggesting a loss of spatial coherence. The building blocks of this analysis 
are the heteroclinic orbits used to construct the typical solitary waves mediating Gunn-like 
oscillations. During most of the oscillation, the motion of the heteroclinic orbits and the 
change of the electric field inside and outside the solitary waves (enclosed by heteroclinic 
orbits) follow adiabatically the evolution of the total current density. When a solitary wave 
reaches the receiving contact, the current increases abruptly and the asymptotic model ade- 
quately approximates this as instantaneous. As an outcome, we have found a criterion that 
shows that single or multiple wave shedding is possible during each oscillation, depending 
on the resistivity of the injecting contact. While single shedding is the usual (stable) Gunn 
effect, multiple wave shedding may break the spatial coherence of the electric field within 
the sample. This new instability mechanism provides an explanation for the complicated 
behavior observed in experiments performed in long semiconductor samples and in the 
numerical simulation of the drift-diffusion model. We have confirmed these results by direct 
numerical simulation of the reduced model, in particular the new predictions of multiple 
shedding of solitary waves in the unstable case. Although simulations and analyses have 
been presented here for the p-Ge model, the general approach is quite general and applies 
to a wide class of pde models which possess the following common properties: 1) an integral 
(over space) constraint; 2) standard boundary conditions which permit multiple stationary 
states - i.e., negative differential resistance; and 3) solitary waves (i.e., pulses) and fronts. It 
is interesting to speculate that different models may lead to the same class of long-sample 
asymptotic limiting model which completely determines and explains the long time dynam- 
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ics (including chaotic temporal behavior and loss of spatial coherence) of the respective full 
models. 
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FIGURES 



FIG. 1. Stationary homogeneous current density, j{E), and linear contact characteristic for a 
sample with the contact resistivity, p =780 f2 cm and compensation ratio Q;=(acceptor concentra- 
tion)/donor concentration) =1.21. 

FIG. 2. Space-time evolution of the electric field E{x, t) and the corresponding current density 
J(t), with parameter values (j) = 6.25 V/cm and p = 780 O cm. The gray scale ranges from 5.3 
V/cm (black) to 14.0 V/cm (white). 

FIG. 3. Space-time evolution of the electric field and the corresponding current for = 7.25 
V/cm. 

FIG. 4. Space-time evolution of the electric field and the corresponding current for a chaotic 
state with = 10.0 V/cm and p = 780 $7 cm. 

FIG. 5. DC bias bifurcation diagram (a) local maxima in the current and (b) largest dimen- 
sionless Lyapunov exponent. 

FIG. 6. This figure shows the velocities of the heteroclinic orbits between E\{J) and E^^i^J), c+ 
and E^{J) and Ei{J), c_ vs. J, both of them in dimensionless units. 

FIG. 7. Space-time evolution of the electric field and the corresponding current density de- 
termined from the asymptotic model for the period-2 state with = 7.25 V/cm and Jc = 9.983 
mA/cm^. 

FIG. 8. Space-time evolution of the electric field and the corresponding current density deter- 
mined from the asymptotic model and showing a chaotic state for parameters Jc = 9.983 mA/cm^ 
and (p = 10.0 V/cm. 
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